Demo guide for package minplascalc

Welcome to minplascalc, a Python 3 package to assist with calculations of equilibrium compositions and thermophysical properties of thermal plasmas of the sort typically encountered in metallurgical processing applications.

Target audience

Plasma technology researchers and professionals with a basic knowledge of the Python programming language.

Module Description

Ionised gases or plasmas are used in many industrial applications such as arc welding, plasma spraying, and electric furnace metallurgy. Engineering plasmas typically operate at atmospheric pressures and temperatures of the order of 104 K. Thermal plasmas of the sort considered here are assumed to be in local thermodynamic equilibirum (LTE), meaning that a single unique temperature can be used to describe them (as opposed to two-temperature plasmas in which the electron temperature is distinct from the heavy-particle temperature). A major advantage of the LTE approximation is that all thermophysical properties of an equilibrium mixture of an arbitrary number of plasma species can be expressed as (complicated) functions of temperature and pressure only - composition is implicit since it is uniquely determined by the state parameters.

Knowledge of these thermophysical properties is of great value to engineers working in plasma technology. Such information is useful for simple design calculations, and is necessary input data for computational fluid dynamics and magnetohydrodynamics models of plasma systems. The calculation of composition and thence the thermophysical properties of a thermal plasma given some fundamental information about the species present is a well-understood but mathematically and numerically complex process. It is prone to error if performed manually, hence the need for this tool.

Things you can calculate with minplascalc: * Statistical mechanics partition functions for individual species using information about the energy levels and excited states * Equilibrium plasma composition in terms of number densities of a specified mixture of species, using the principle of Gibbs free energy minimisation at a specified temperature and pressure * Density \(\rho\), (relative) enthalpy \(H\), and heat capacity \(C_P\) at equilibrium conditions

Things you can’t calculate with minplascalc: * Compositions or thermophysical properties of two-temperature plasmas * \(C_P\) for arbitrary non-LTE compositions * Plasmas of complex molecules or molecular ions consisting of more than two atoms * Plasmas involving negatively-charged ions * Transport or thermal radiation properties (see future versions though)

The package also includes some basic tools to format and import data for individual species obtained from the NIST Atomic Spectra and Chemistry WebBook databases, and store it for use by future simulations.

Partition functions

The starting point for thermal plasma calculations is generally the statistical mechanics partition functions for each of the species present. Users of minplascalc should not normally need to access these functions explicitly as they are incorporated directly into the composition and thermophysical property calculators, but they are exposed in the API in case the need to do so ever arises.

Recall that the partition function for a particular species is a description of the statistical properties of a collection of atoms or molecules of that species at thermodynamic equilibrium. Partition functions are normally presented as the sum of weighted state probabilities across the possible energy states of the system. In general at moderate plasma temperatures up to a few 104 K, a species’ total partition function \(Q_{tot}\) can be written as the product of several unique partition functions arising from different quantum mechanical phenomena (assuming weak state coupling and no contribution from nuclear states):

\(Q_{tot} = Q_t Q_{int} = Q_t Q_e Q_v Q_r\)

Here, \(Q_t\) is the translational partition function due to the species’ ability to move around in space, \(Q_{int}\) is the internal partition function due to energy states internal to the particles of the species, \(Q_e\) is the electronic partition function due to different possible arrangements of the electron structure of the species, \(Q_v\) is the vibrational partition function due to the ability of the bonds in a polyatomic species to vibrate at different energy levels, and \(Q_r\) is the rotational partition function due to a species’ ability to rotate around its center of mass at different energy levels.

minplascalc distinguishes between three different types of species - monatomic (charged or uncharged single atoms), diatomic (charged or uncharged bonded pairs of atoms), and free electrons. The formulae used for the various partition functions for each are shown in the table below.

Partition Function

Monatomic

Diatomic

Electron

\(Q_t\), m-3

\[{\left ( \frac{2 \pi m_s k_B T}{h^2}\right )}^{\frac{3}{2}}\]
\[{\left ( \frac{2 \pi m_s k_B T}{h^2}\right )}^{\frac{3}{2}}\]
\[{\left ( \frac{2 \pi m_e k_B T}{h^2}\right )}^{\frac{3}{2}}\]

\(Q_e\), dim’less

\[\sum_j g_j \exp \left(-\frac{E_j}{k_B T}\right)\]
\[g_0\]
\[2\]

\(Q_v\), dim’less

\[1\]
\[\frac{1}{1-\exp\left( -\frac{\omega_e}{k_B T} \right)}\]
\[1\]

\(Q_r\), dim’less

\[1\]
\[\frac{k_B T}{\sigma_s B_e}\]
\[1\]

Here \(m_s\) and \(m_e\) are the mass of one particle of the species concerned, \(k_B\) is Boltzmann’s constant, \(T\) is temperature, \(h\) is Planck’s constant, \(g_j\) and \(E_j\) are the quantum degeneracy and energy (in J) of electronic energy level j (with j = 0 being the ground state), and \(\omega_e\), \(\sigma_s\) and \(B_e\) are the vibrational, symmetry, and rotational constants respectively for a diatomic molecule.

Ionisation energy lowering

In general the ionisation energy required to remove a single electron from a particle of a species is a constant for that particular species when considered in isolation. However, in a mixture of different species and free electrons, the ionisation energy is lowered by a small amount due to local electrostatic shielding effects. This affects both the calculation of the partition functions (the summation of electronic state contributions for monatomic species ignores states with energies above the lowered ionisation energy) and the calculation of equilibrium plasma compositions (the equilibrium relationships are defined using the reference energy levels for each species, which in turn depend on the lowered ionisation energies). Ionisation energy lowering is a complex problem in plasma physics, but there exist many approximate methods for quantifying this effect using the theory of Debye-shielded potentials. Provided the same method is used for all species, the calculation errors generally remain small. The ionisation energy lowering calculation is not exposed to the user in the minplascalc API, since it is only required internally for calculation of species partition functions and LTE compositions.

minplascalc uses the analytical solution of Stewart and Pyatt 1966 (see references in README). In this method, the ionisation energy lowering for each species is calculated explicitly using:

\[\frac{\delta E_i}{k_B T} = \frac{\left [ \left (\frac{a_i}{\lambda_D} \right )^3 + 1 \right ]^\frac{2}{3} -1}{2 \left( z^*+1 \right)}\]

where:

\[z^* = \left ( \frac{\sum z_j^2 n_j}{\sum z_j n_j} \right )_{j \neq e}\]
\[a_i = \left ( \frac{3 \left (z_i + 1 \right )}{4 \pi n_e} \right )^\frac{1}{3}\]
\[\lambda_D = \left ( \frac{k_B T}{4 \pi e^2 \left ( z^* + 1 \right ) n_e} \right )^\frac{1}{2}\]

Here, \(\delta E_i\) is the amount the ionisation energy of species i is lowered by (in J), \(a_i\) is the ion-sphere radius of species i, \(\lambda_D\) is the Debye sphere radius, \(z^*\) is the effective charge number in a plasma consisting of a mixture of species of different charges, \(z_j\) is the charge number of species j, \(n_j\) is the number density (particles per cubic meter) of species j, and \(e\) is the electron charge.

Calculation of LTE compositions

Given temperature, pressure, and a set of species present in a plasma (and some information about the elemental composition of the mixture if more than one element is present), the number density of each species at thermodynamic equilibrium can be calculated using the principle of Gibbs free energy minimisation. This is an important intermediate step in calculating the thermopysical properties, and may also be useful in its own right if one is interested in the relative proportions of different species in complex plasmas. It is exposed to the user in the minplascalc API.

To start, recall the definition of Gibbs free energy:

\[G = G^0 + \sum_i \mu_i N_i\]

where \(G\) is the Gibbs free energy of a system, \(G^0\) is a reference value depending only on temperature and pressure, \(\mu_i\) is the chemical potential of species i, and \(N_i\) is the absolute number of particles of species i present. In terms of statistical mechanics properties, \(\mu_i\) can be represented as:

\[\mu_i = E_i^0 - k_B T \ln \left ( \frac{Q_{tot,i}V}{N_i} \right )\]

where \(Q\) is the partition function defined earlier, \(E_i^0\) is the reference energy of the species relative to its constituent uncharged atoms (for uncharged monatomic species and electrons \(E_i^0=0\), for uncharged diatomic species it is the negative of the dissociation energy, and for charged species it is \(E_i^0\) of the species with one fewer charge number plus the lowered ionisation energy of that species), and \(V\) is the volume of the system. From the ideal gas law, we have:

\[V = \frac{k_B T \sum_i N_i}{p}\]

where \(p\) is the specified pressure of the system.

A system at equilibrium is characterised by a minimum stationary point in \(G\), giving an independent equation for each species i which simplifies to:

\[\frac{\partial G}{\partial N_i} = \mu_i = 0\]

This set of equations must be solved subject to constraints supplied by the conservation of mass of each element present:

\[\sum_i v_{ij} N_i = \eta_j^0\]

where \(v_{ij}\) is the stoichiometric coefficient representing the number of atoms of element j present in species i, and \(\eta_j^0\) is the (fixed) total number of atoms of element j present in the system, obtained from user specifications. Together with this, one additional constraint is supplied by the requirement for electroneutrality of the plasma:

\[\sum_i z_i N_i = 0\]

In minplascalc, the previous three sets of equations are solved using an iterative Lagrange multiplier approach to obtain the set of \(N_i\) (and hence number density \(n_i = N_i / V\)) at LTE starting from an initial guess.

Calculation of plasma density

Given a plasma composition in terms of number densities \(n_i\), the mass density is a straightforward calculation:

\[\rho = \frac{1}{N_A} \sum_i n_i M_i\]

where \(M_i\) is the molar mass of species i in kg/mol, and \(N_A\) is Avogadro’s constant. The density calculation is exposed as a function call in the minplascalc API.

Calculation of plasma enthalpy

Calculation of the plasma enthalpy at a particular temperature, pressure, and species composition is performed using the statistical mechanics definition of internal energy:

\[U = -\sum_j \frac{1}{Q_j} \frac{\partial Q_j}{\partial \tau}\]

where \(U\) is the internal energy in J/particle for a particular species, \(Q_j\) are the various kinds of partition functions making up \(Q_{tot}\) for the species, and \(\tau=(k_B T)^{-1}\). Formulae for \(U\) of various plasma species are thus readily produced using the expressions for \(Q_j\) given earlier.

Recall the thermodynamic definition of enthalpy:

\[H = U + p V\]

When multiple species are present, the relative reference energy \(E_i^0\) for each species must also be included. Application of the ideal gas law to the \(pV\) term then gives:

\[H_i = U_i + E_i^0 + k_B T\]

where \(H_i\) is the enthalpy of species i in J/particle. Summing over all component species of a plasma and dividing by the density then gives the total enthalpy of the mixture in J/kg:

\[H = \frac{\sum_i n_i H_i}{\rho} = N_A \frac{\sum_i n_i H_i}{ \sum_i n_i M_i}\]

The enthalpy calculation is exposed to the user in the minplascalc API via a function call, however, it is important to note that in the current formulation some values of \(E_i^0\) may be negative and the calculated enthalpy value will therefore be relative to an arbitrary non-zero value.

Calculation of plasma heat capacity

A direct calculation of \(C_P\) given an arbitrary, non-LTE plasma composition is possible if some knowledge of the reaction paths between species is also supplied. Although any set of consistent reaction paths will give the same result, choosing one actual set of paths from the many possible options implies that it represents reality, and this is certainly open to some debate. In the spirit of keeping minplascalc focused on path-independent equilibrium plasma problems, the heat capacity calculation is instead performed by numerical derivative of the enthalpy around the temperature of interest:

\[C_P = \left( \frac{\partial H}{\partial T} \right)_p = \frac{H_{T+\Delta T,p} - H_{T-\Delta T,p}}{2 \Delta T}\]

Here, \(H_{T+\Delta T,p}\) and \(H_{T-\Delta T,p}\) are enthalpy calculations for the LTE plasma composition at fixed pressure, and temperatures slightly above and slightly below the target temperature \(T\). This calculation is exposed to the user in the minplascalc API via a function call, and it is important to note that it only gives the heat capacity of LTE compositions.